# Seurat official vignette
# https://satijalab.org/seurat/articles/visualization_vignette.html

library(Seurat)
library(tidyverse)
library(ggpubr)
data("pbmc3k.final")
pbmc3k.final$groups <- sample(c("group1", "group2"), size = ncol(pbmc3k.final), replace = TRUE)
features <- c("LYZ", "CCL5", "IL32", "PTPRCAP", "FCGR3A", "PF4")
pbmc3k.final

RidgePlot(subclstr10x[["hB01_PlasmaB-IgG"]], features = 'FCGR2B')
VlnPlot(subclstr10x[["hB01_PlasmaB-IgG"]], features = 'FCGR2B')

select10x <- subclstr10x[c("hB01_PlasmaB-IgG",
                           "hB04_FollicularB-MS4A1",
                           'hM03_cDC2-CD1C',
                           'hM09_Macro-PLTP',
                           'hM10_Macro-IL1B',
                           'hM11_Monolike-FCN1',
                           'hM12_TAM-C1QC',
                           'hM13_TAM-SPP1')]
slct10x <- reduce(select10x, merge)
VlnPlot(labelSmart, group.by = 'Sub_Cluster',
        features = 'FCGR2B') + theme(legend.position = "none")
VlnPlot(label10x,
        features = 'FCGR2B', 
        group.by = 'Sub_Cluster',
        split.by = 'ITgeno',
        split.plot = TRUE)

VlnPlot(labelSmart,
        features = 'FCGR2B', 
        group.by = 'Sub_Cluster',
        split.by = 'ITgeno',
        split.plot = TRUE)

# Visualize the number of cell counts per sample
label10x@meta.data %>% 
  ggplot(aes(x=Sub_Cluster, fill = ITgeno)) + 
  geom_bar() +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  theme(plot.title = element_text(hjust=0.5, face="bold")) +
  ggtitle("NCells")

DotPlot(label10x,
        features = 'FCGR2B', 
        #group.by = 'Sub_Cluster',
        split.by = 'ITgeno')
# Single cell heatmap of feature expression
DoHeatmap(subset(label10x, downsample = 100),
          features = 'FCGR2B',
          group.by = 'Sub_Cluster',
          slot = 'data',
          size = 3) + theme(legend.position = "none")


# KEGG pathway hierarchy ---------
library(KEGGREST)

gel.pbs.kegg <- read_csv('/home/Bill/Downloads/Gel_vs_PBS.KEGGSummary.csv')

kegg.ip.subset <- read_delim('Archive/covid19/ref/kegg_info_process.txt', col_names = c('id','name'))

kegg.immu.subset <- read_delim('Archive/covid19/ref/kegg_immune.txt', col_names = c('id','name'))

kegg.lcf.subset <- bind_rows(kegg.immu.subset, kegg.ip.subset, .id = 'subclass') |>
  mutate(subclass = case_match(subclass,
                               '1' ~ 'Immune System',
                               '2' ~ 'Environmental Information Processing'))

gel.pbs.ip <- gel.pbs.kegg |>
  mutate(id = str_remove(PathwayID, 'mmu') |> as.numeric()) |>
  right_join(kegg.lcf.subset) |>
  filter(!is.na(PathwayID) & qvalue < .05) |>
  separate(BgRatio, into = c('bg_set', 'bg_all'), convert = TRUE) |>
  mutate(rich_factor = Count / bg_set) |>
  write_csv(str_glue('{dir_download}lcf_kegg_subset.csv'))

gel.pbs.ip |>
  filter(str_detect(subclass, 'Immune')) |>
  slice_max(rich_factor, n = 20) |>
  ggplot(aes(fct_reorder(Description, rich_factor), rich_factor, color = -log10(qvalue), size = Count)) +
  geom_point() +
  geom_text(aes(label = Count), color = 'black', size = 4) +
  labs_pubr() +
  theme_pubr(legend = 'right') +
  coord_flip() +
  scale_color_gradient(low = 'yellow', high = 'red') +
  scale_size_area(max_size = 10) +
  labs(title = str_wrap('20 most enriched KEGG pathways in "Immmune System" subclasses', width = 40),
       x = 'Pathway Description',
       color = '-log10(FDR)')

ggsave(path = dir_download, 'lcf_immune_kegg.pdf', width = 10, height = 7.5)

gel.pbs.ip |>
  filter(str_detect(subclass, 'Environmental')) |>
  slice_max(rich_factor, n = 20) |>
  ggplot(aes(fct_reorder(Description, rich_factor), rich_factor, color = -log10(qvalue), size = Count)) +
  geom_point() +
  geom_text(aes(label = Count), color = 'black', size = 4) +
  labs_pubr() +
  theme_pubr(legend = 'right') +
  coord_flip() +
  scale_color_gradient(low = 'yellow', high = 'red') +
  scale_size_area(max_size = 10) +
  labs(title = str_wrap('20 most enriched KEGG pathways in "Environmental Information Processing" subclasses', width = 40),
       x = 'Pathway Description',
       color = '-log10(FDR)')

ggsave(path = dir_download, 'lcf_signal_kegg.pdf', width = 10, height = 7.5)
